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Abstract.  This  paper  presents  an  approach  to  modeling  discrete  features  on  a  U.S.  Army 
endurance  test  course.  The  features  in  this  case  are  39  moguls  with  heights  varying  between 
0.36  and  0.63  m  built  into  a  portion  of  the  course.  Because  they  are  enumerable  and  localized , 
they  do  not  lend  themselves  to  traditional  modeling  techniques  such  as  RMS,  PSD,  IRI,  etc.  We 
assume  that  the  mogul  has  a  well  defined  shape  and  we  evaluate  five  different  approaches  to 
modeling  this  shape  assuming  that  each  mogul  may  be  thought  of  as  an  ideal  shape 
superimposed  with  stationary  noise.  The  models  that  we  evaluate  are  the  Gaussian  function, 
the  Hanning  window,  the  support  vector  machine,  a  two-Gaussian  function  and  a  two-Gaussian 
function  with  fixed  width  and  amplitude  ratios.  We  evaluate  these  models  by  computing  the 
RMS  error  of  each  best  fit  and  by  evaluating  the  stationarity  of  the  residue.  We  further  evaluate 
each  mogul  model  by  running  a  HMMWV  dynamics  model  over  each  and  comparing  several 
responses  to  those  obtained  from  the  profiled  mogul.  In  our  analysis  we  find  that  a  model 
consisting  of  two  Gaussian  functions  with  related  widths  and  amplitudes  yields  an  unbiased 
estimate  of  a  mogul  and  can  be  made  to  approximate  any  mogul  by  adjusting  its  width  and 
amplitude. 
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HMMWV. 


1  Introduction 

The  U.S.  Army  performs  endurance  testing  of  its  vehicles  at  several  proving  grounds  each  of 
which  have  various  courses  to  represent  primary,  secondary  and  cross-country  terrains.  These 
courses  are  used  by  the  Army  to  evaluate  vehicles  against  contractual  reliability  requirements.  It 
is,  therefore,  in  the  interest  of  both  the  Army  and  its  vehicle  developers  to  measure,  characterize, 
and  understand  the  essential  features  of  these  courses.  Such  an  understanding  helps  the  Army 
monitor  and  maintain  course  severity  and  helps  the  contractor  know  how  its  vehicle  will  be 
evaluated.  All  terrain  monitoring  begins  with  a  measurement  of  the  terrain  profile,  for  which 
there  are  many  available  techniques  which  are  nicely  summarized  by  Sayers  and  Karamihas 
(1998).  A  road  profile  is  a  measurement  of  terrain  elevation  vs.  horizontal  distance  along  the 
course.  The  horizontal  resolution  of  such  profiles  varies  from  0.3  m  to  1  mm;  they  may  have 
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lengths  varying  from  30  m  to  10  km  (or  more)  and  are  typically  measured  for  two  vehicle  tracks. 
The  volume  of  information  and  the  random  nature  of  the  raw  profile  data  motivates  their 
reduction  to  some  essential  features  which  are  preferably  invariant  to  the  statistical  variation 
observed  in  the  profile  measurement.  It  is  desired  that  these  features  correlate  to  a  vehicle 
performance  metric  such  as  ride  quality  or  durability.  A  one  statistic  traditionally  used  by  the 
Army  has  been  the  root  mean  square  (RMS)  of  the  detrended  terrain  profile  which  works  well  for 
terrains  which  are  stationary.  It  is  appealing  because  it  captures  the  roughness  in  one  number. 
Recently,  the  Amy  has  begun  to  characterize  its  profiles  with  the  power  spectrum  density  (PSD) 
(or  more  correctly  the  wave  number  spectrum).  Typically  these  PSDs  are  then  used  to 
generalize  the  frequency-domain  shape  to  a  sloped  line  which  approximates  the  PSD.  Such 
representations  are  usually  documented  in  the  form  of  a  frequency-domain  plot  with  labeled 
slopes  and  breakpoints.  The  PSD  is  again  appealing  in  that  it  captures  the  essential  terrain 
information  in  a  few  parameters  represented  in  a  plot;  it  however  assumes  that  the  terrain  profile 
is  statistically  stationary  (see  Steinwolf  et  al,  2005).  The  automotive  industry  has  developed 
their  own  metrics  for  road  roughness  called  the  International  Roughness  Index  (IRI),  and  the 
Speed  Roughness  Index  (SRI).  The  IRI  is  a  computation  based  on  the  cumulative  suspension 
displacement  of  the  so  called  “Golden  car  model”  run  over  the  profile  at  a  fixed  speed.  The  SRI 
adjusts  the  IRI  for  speed. 

Given  these  characterizing  features  (e.g.  RMS,  PSD,  IRI),  it  is  often  desired  to  reverse  this 
process;  that  is,  given  some  roughness  metric,  synthesize  a  terrain  profile  which  is  a  faithful 
representation  of  the  original  terrain  profile.  In  the  synthesis  problem,  it  is  important  that  enough 
features  be  retained  from  the  original  profile  to  allow  proper  reconstruction.  Single-value 
metrics  such  as  RMS,  IRI,  and  SRI  are  not  suited  to  profile  synthesis.  The  PSD  feature  of  a 
profile  has  been  used  by  the  Army  in  terrain  synthesis,  however,  the  resulting  profiles  are 
stationary  and  have  random  phase  relationships,  normally  not  like  the  actual  profile.  Sandu  et 
al  (2005)  has  proposed  terrain  synthesis  based  on  polynomial  chaos.  Sandu  and  Li  (2007) 
recently  proposed  terrain  synthesis  based  on  ARMA  and  ARIMA  models.  Sun  et  al  (2005) 
applied  uniform  modulated  processes  (UMP)  and  empirical  mode  decomposition  (EMD)  to  the 
synthesis  of  non-stationary  terrain  profiles  such  as  are  typically  measured  at  Army  proving 
grounds. 

Each  of  these  methods  of  terrain  characterization  and  synthesis  approach  the  problem  in  an 
agnostic  manner,  treating  the  terrain  profile  as  a  time-series  with  certain  assumed  properties  (i.e. 
Gaussian,  stationary,  correlated,  etc).  None  of  them  however,  address  the  circumstance  in 
which  there  is  a  mix  of  random  and  non-random  features  in  a  terrain  profile.  Although  this  is 
normally  a  valid  assumption,  it  is  not  always  the  case.  The  topic  of  this  paper  is  to  present  an 
approach  to  modeling  particular  non-random  features  (moguls)  on  a  particular  Army  proving 
ground  course  for  the  purposes  of  both  characterization  and  synthesis.  It  is  the  intent  of  this 
paper  to  statistically  characterize  these  moguls  using  a  number  of  candidate  functions  and 
determine  which  is  most  suitable.  We  choose  direct  metrics  based  on  quality  of  fit  and  indirect 
metrics  derived  from  HMMWV  dynamics  model  traversing  the  moguls.  Functions  which  are  fit 
to  the  moguls  are  Gaussian,  Flanning,  Support  Vector  Machines,  two  summed  Gaussians,  and 
to  ratio-related  Gaussians. 

2  Profile  Description 

The  particular  course  in  view  is  used  by  the  Army  for  endurance  testing  and  has  a  particular 
section  containing  39  moguls  with  heights  varying  between  0.35  and  0.65  m  spaced  at 
approximate  30  m  intervals.  Figure  1  shows  a  section  of  this  terrain  which  attempts  to  place  the 
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Figure  1.  Three  plots  of  the  course  being  analyzed,  (left)  2  km  section  of  the  course 
containing  the  39  moguls,  (middle)  view  of  a  single  pass  of  moguls,  (right)  view  of  three 
moguls  on  a  representative  grade.  The  plots  are  not  to  scale. 


moguls  in  reference  to  the  entire  profile.  The  profile  data  used  in  the  following  analysis  was 
taken  from  this  course  in  2002.  The  profile  consists  of  correlated  right  and  left  track  separated 
by  1.22  m  (4  ft)  and  the  longitudinal  resolution  is  7.62cm  (3  in.).  The  device  which  measured 
the  profile  was  developed  by  the  U.S.  Army  Aberdeen  Proving  Ground  (APG)  for  the  monitoring 
of  course  roughness  (Castaldo  et  al,  1991). 

3  Modeling  Approaches 

To  begin  we  note  that  most  of  the  moguls  are  situated  on  graded  terrain  and  since  we  are 
interested  in  modeling  the  ride  and  shock  aspects  of  them  we  will  remove  the  grade  information 
in  our  analysis.  We  selected  each  mogul  by  finding  its  center  point.  Each  mogul  was  extracted 
from  the  terrain  profile  by  taking  120  sample  points  ahead  and  behind  this  center  point,  so  each 
mogul  was  roughly  contained  at  the  center  of  a  18.28  m  (60  ft)  stretch  of  terrain  (segment). 
Each  segment  was  then  subjected  to  a  linear  least  squares  analysis  and  the  prevailing  bias  and 
slope  were  removed.  A  plot  of  these  detrended  mogul  segments  is  shown  in  Figure  2.  Next  we 
chose  the  point  in  the  horizontal  axis  where  we  regard  the  mogul  to  begin  and  end.  By 
observing  Figure  2  it  is  clear  that  the  traces  of  the  moguls  agree  fairly  well  in  the  interval  [-6m, 
+6m]  and  that  the  agreement  between  the  moguls  lessens  at  distances  more  than  6  m  from  the 
center.  We  will  therefore  perform  our  analysis  on  an  interval  of  ±  6m  (±  20  ft)  resulting  in  an 
assumed  width  of  12.2  m  (40  ft).  It  is  clear  in  Figure  2  that  the  heights  of  the  moguls  vary  from 
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Figure  2.  Plot  of  39  moguls  overlaid,  (left)  plot  over  the  segment  [-9m,  9m],  (right)  plot 
over  the  horizontal  interval  [-6m,  6m]  for  comparison  to  the  fit  moguls.  Right  track  traces 
are  in  blue  and  left  track  traces  are  in  green.  Observe  that  the  moguls  range  in  height 
from  0.36  m  to  0.63  m. 
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approximately  0.36  m  to  0.63m  and  we  expect  our  analysis  to  show  that  the  widths  vary  as  well. 

In  our  modeling  approaches  we  will  seek  to  minimize  the  mean  squared  error  (MSE)  between 
the  actual  mogul  data  and  our  model.  We  will  regard  each  mogul  as  a  set  of  ordered  pairs 

{(X/’37I)}|  where  xt  is  the  horizontal  distance,  y.  is  the  corresponding  elevation  and  i  is  the 

number  of  samples,  which,  in  our  case  is  161.  Our  models  will  take  the  continuous  form 
y  =  f(x; p)  where  x  is  the  independent  variable,  p  is  a  vector  of  parameters  to  be  optimized. 
We  use  the  “hat”  notation  to  denote  an  estimated  value.  Finally  we  define  the  particular 
estimates  as  yi=f(xi;p)  .  Our  optimization  problem  will  then  consist  of  minimizing  the 
functional 


min:  Z(/(x,;p)-y,T  +  cp(p)  (i) 

where  V(-)  is  a  regularization  functional  and  C  is  a  regularization  constant  (which  we  will 

normally  choose  to  be  zero).  Then  by  observing  Figure  2  we  are  inclined  be  believe  that  the 
shape  of  these  moguls  closely  resembles  a  bell  type  curve,  so  our  first  attempts  will  be  based 
on  the  Gaussian  function  and  the  Flanning  window. 

3. 1  The  single  Gaussian  function 

We  will  begin  our  search  for  the  best  fit  of  the  mogul  by  choosing  a  Gaussian  function  of  the 
following  form 

-(x-x? 

f(x\p)  =  ae  2ff2  +bx  +  c  (2) 

where  the  vector  of  free  parameters  is  p  =  [«  x  a  b  r]7  .  Although  the  position  x,  bias  c 

and  slope  b  are  products  of  the  optimization,  they  are  present  to  eliminate  biased  errors  due  to 
shifts  in  the  reference  data.  In  this  representation,  we  will  take  the  value  a  to  be  the  mogul 
amplitude  and  since  <x  is  the  standard  deviation,  we  will  regard  the  width  of  the  bump  as 
w  =  6(7  (i.e.  ±3cr ).  After  optimizing  (1 )  given  (2)  the  fit  functions  are  shown  in  Figure  3.  There  it 
may  be  observed  that  the  error  is  biased  based  on  the  distance  from  the  center  of  the  mogul. 
The  RMS  errors  for  the  right  track  are  on  the  interval  [0.95  cm,  2.58  cm]  with  an  average  of  1 .63 
cm.  Likewise,  for  the  left  track  the  RMS  error  falls  in  the  range  [0.79  cm,  2.27  cm]  with  an 
average  of  1 .62  cm.  This  bias  in  the  error  motivates  us  to  look  for  a  better  candidate  function. 
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Figure  3.  (left)  Plot  of  the  Gaussian  best  fit  for  the  39  moguls,  (middle)  Plot  of  the  mogul 
amplitude  {a)  vs.  the  mogul  width  (6a).  (right)  Plot  of  the  error  from  the  Gaussian  best 
fit.  Note  scales  are  not  the  same. 
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3.2  The  Hanning  window  function 


Our  next  candidate  function  is  the  so  called  Hanning  window  which  has  been  used  on  some 
occasions  as  a  road  disturbance  input  (see  Lin  et  al,  1995).  In  our  case  the  Hanning  function 
takes  the  following  form 


/(*;p) 


0, 


otherwise 


(3) 


where  the  vector  of  free  parameters  is  p  =  [a  x  w  b  c\  .  In  this  representation,  we  will  take 
the  value  a  to  be  the  mogul  amplitude  and  w  as  the  width.  After  optimizing  (1)  given  (3)  the  fit 
functions  are  shown  in  Figure  4.  There  it  may  be  observed  the  error  is  biased  based  on  the 
distance  from  the  center  of  the  mogul  as  it  was  for  the  Gaussian  function.  The  RMS  errors  for 
the  right  track  are  on  the  interval  [1.29  cm,  3.18  cm]  with  an  average  of  2.33  cm.  Likewise,  for 
the  left  track  the  RMS  error  falls  in  the  range  [1 .33  cm,  3.08  cm]  with  an  average  of  2.28  cm. 
Examining  these  results  with  those  of  the  Gaussian  best  fit  it  appears  that,  in  general,  the 
Gaussian  function  is  a  better  representation  of  the  moguls  than  the  Hanning  window,  we 
therefore  will  eliminate  the  Hanning  function  from  further  consideration. 


3.3  The  Support  Vector  Machine  (SVM)  function 

Since  our  analysis  thus  far  has  shown  that  one  Gaussian  function  provides  a  reasonable  fit,  we 
speculate  that  a  set  of  superimposed  Gaussian  functions  may  perform  better.  We  therefore  will 
employ  a  recently  developed  technology  called  Support  Vector  Machines  (SVM)  (see  Crisianini 
and  Shawe-Taylor  (2000)  for  a  full  development  of  the  theory)  to  find  a  suitable  class  of 
functions  with  which  to  fit  our  moguls.  The  form  of  the  SVM  estimator  is  as  follows 

/(*;  v)  =  YJPie~y(x~x,)  +b  (4) 

i= 1 

where  the  vector  of  free  parameters  is  p  =  [/?,  ■■■  j3f  bf  consisting  of  the  weights  /?. 

and  the  bias  b  and  the  values  x.  are  the  horizontal  distances  from  our  mogul  center.  In  SVM 

theory  k(x,y)  =  e~r(-x~y)  is  called  the  Gaussian  kernel  function  for  which  we  are  free  to  choose 
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Figure  4.  (left)  Plot  of  the  Hanning  window  best  fit  for  the  39  moguls,  (middle)  Plot  of 
the  mogul  amplitude  (a)  vs.  the  mogul  width  (w).  (right)  Plot  of  the  error  from  the 
Hanning  window  best  fit.  Note  scales  are  not  the  same. 
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the  kernel  parameter  y .  The  parameters  p  are  found  by  solving  the  following  quadratic 
programming  problem 


min : 

(A) 


subject  to: 


ZZAA/’"-"’ 


o  ' 

i=l  7=1  1=1 

ZA=o,  -c<p<c 


& 


(5) 


where  e  indicates  an  upper  absolute  bound  on  the  error  we  wish  to  tolerate  and  C  is  a 
regularization  constant.  In  our  analysis  we  chose  ^=0.149  1/m2  (0.0138  1/ft2)  as  the  kernel 
parameter,  which  corresponds  to  a  width  of  10.97  m  (36  ft)  and  was  determined  by  examining 
the  average  width  generated  in  the  Gaussian  fit  (Figure  3).  We  chose  the  error  bound  e  to  be 
2.54  cm  (1  inch)  and  the  regularization  constant  to  be  63.5  cm  (25  inches).  By  design,  the 

solution  of  (5)  is  sparse  in  the  set  of  weights  {#}|  and  those  x(  for  which  ft* 0  are  called 
support  vectors.  Also  the  value  of  the  bias  b  is  uniquely  determined  by  the  so  called  KKT 
conditions  and  the  values  of  .  In  the  optimization  it  was  found  that  the  number  of  support 

vectors  varied  from  5  to  16.  After  optimization  of  (5)  the  fit  functions  are  shown  in  Figure  5.  We 
observe  that  the  error  for  the  most  part  is  unbiased  and  is  mostly  clamped  at  ±e,  which  the 
SVM  algorithm  does  by  design.  The  RMS  errors  for  the  right  track  are  on  the  interval  [1.27  cm, 
1 .63  cm]  with  an  average  of  1 .42  cm.  Likewise,  for  the  left  track  the  RMS  error  falls  in  the  range 
[1.20  cm,  1.65  cm]  with  an  average  of  1.42  cm.  Comparing  these  results  to  both  the  Gaussian 
Flanning  window  approaches,  the  SVM  clearly  out  performs  them  both  in  terms  of  average  RMS 
error  however,  the  Gaussian  function  does  better  for  some  of  the  moguls  as  indicated  by  the 
minimum  RMS  error.  The  strength  of  the  SVM  approach  is  that  by  adjusting  the  parameter  e 
we  can  obtain  an  (almost)  arbitrary  amount  of  precision.  The  weakness  of  the  SVM  approach  is 
that  it  is  a  black  box  algorithm  that  yields  no  particular  insight  to  the  generalized  shape  of  a 
typical  mogul.  Since  our  goal  is  find  a  mogul  model  with  two  parameters,  amplitude  and  width, 
we  will  continue  to  evaluate  other  models. 


3.4  The  Two-Gaussian  function 

Knowing  that  the  superposition  of  several  Gaussians  improves  the  fit  (at  least  for  some  moguls) 
we  will  next  build  a  hypothesis  function  composed  of  two  Gaussian  functions.  These  functions 


-5  0  5 

Horizontal  Distance  [m] 


E_ 

o 

% 

g5T 


-5  0  5 

Horizontal  Distance  [m] 


0.1 

0.05 


Horizontal  Distance  [m] 


Figure  5.  (left)  Plot  of  the  SVM  best  fit  for  the  39  moguls,  (middle)  Plot  of  the  weights  /3i 
vs.  the  support  vectors  xi  for  /? .  *0.  (right)  Plot  of  the  error  from  the  SVM  best  fit.  Note 
scales  are  not  the  same. 
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will  have  freedom  to  vary  independently  in  amplitude  and  width,  however  they  will  be  centered 
on  the  same  point,  x .  Our  candidate  function  takes  the  following  form 

-(x-xf  -{x-xf 

f(x',p)=ale  2cr‘  +a2e  2cTl  +bx+c  (6) 

where  the  vector  of  free  parameters  is  p  =  a2  x  at  a2  b  cf .  (It  is  worth  noting  that 

(6)  may  be  thought  of  as  a  1-2-1  artificial  neural  network  with  Gaussian  activation  functions  in 
the  hidden  layer  and  a  linear  output  function.)  In  this  representation,  we  will  take  the  sum 
a1+a2  to  be  the  mogul  amplitude  and  width  we  will  regard  as  w  =  6(alal  +  a2a2)/(al  +  ai)  which 

is  a  weighted  average  of  the  two  Gaussian  widths.  In  this  case,  it  was  discovered  that  the 
optimization  problem  was  ill-posed  (see  Tikhonov  &  Arsenin,  1977)  and  we  therefore  introduced 
the  regularization  functional  V(p)  =  a;  +a]  +  <7,2  +  cr2  and  the  regularization  constant  C  =  0.1 

from  (1)  to  stabilize  the  solution.  After  optimizing  (1)  given  (6)  the  fit  functions  are  shown  in 
Figure  6.  There  it  may  be  observed  that  the  error  is  for  the  most  part  unbiased.  The  RMS 
errors  for  the  right  track  are  on  the  interval  [0.57  cm,  2.07  cm]  with  an  average  of  1.13  cm. 
Likewise,  for  the  left  track  the  RMS  error  falls  in  the  range  [0.58  cm,  2.03  cm]  with  an  average  of 
1.06  cm.  Examining  these  results  with  those  of  the  prior  models  it  appears  that  the 
superposition  of  two  independent  Gaussian  functions  provides  the  best  fit  thus  far.  It  is 
encouraging  to  see  that  the  addition  of  the  second  function  is  able  to  effectively  compensate  for 
the  positional  bias  found  in  the  single  Gaussian  approach.  In  this  case  we  currently  have 
doubled  the  number  of  significant  parameters  from  two  to  four.  We  may  now  ask  if  the 
parameters  a, ,  a2,  cr, ,  and  a2  are  truly  independent.  Shown  in  Figure  7  are  two  plots  which 
will  help  answer  this  question.  The  left  plot  shows  the  width  and  amplitudes  of  both  Gaussians 
independently  which  after  training  the  wider/taller  Gaussians  were  labeled  with  ax  and  cfx  while 
the  narrower/shorter  ones  were  labeled  with  a2  and  a2.  In  this  plot  we  observe  that  the  first 
and  second  Gaussian  parameter  sets  tend  to  cluster  into  two  groups  indicating  that  there  may 
be  a  relationship  between  the  parameters  which  we  will  hypothesize  to  be  proportional.  To  test 
this  hypothesis  we  compute  the  ratios  of  widths  (i.e.  o2/ct{  )  and  amplitudes  (i.e.  a2/al)  for  all 
39  moguls.  These  ratios  are  plotted  in  Figure  7  on  the  right.  Furthermore  we  compute  the 
average  value  of  each  of  these  ratios  to  obtain  the  desired  proportional  relationships  as 
Pa  =  <*2 /°i  =0-431  and  pa  =  a2jax  =0.551 .  Our  next  and  last  set  of  functions  to  evaluate  will 


Figure  6.  (left)  Plot  of  the  two-Gaussian  best  fit  for  the  39  moguls,  (middle)  Plot  of  the 
mogul  amplitude  vs.  width  as  estimated,  (right)  Plot  of  the  error  from  the  two-Gaussian 
best  fit.  Note  scales  are  not  the  same. 
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0.431 


Figure  7.  Plot  of  the  parameter  relationships  for  the  two-Gaussian  best  fit  for  the  39 
moguls,  (left)  The  raw  parameters.  The  a  symbol  represents  6 ax  vs.  a,  (red  =  right, 

black  =  left)  and  the  v  symbol  represents  6<x2  vs.  a2  (blue=right,  green=left).  (right)  The 

parameter  ratios  a2/ax  vs.  cr2/cr,  .  The  red  lines  represent  the  average  value  for  each 
of  these  ratios. 

consist  of  two  Gaussians  with  widths  and  amplitudes  related  by  these  ratios. 


3.5  The  Fixed-Ratio  Two-Gaussian  Function 


Knowing  that  the  superposition  of  two  independent  Gaussians  improves  the  fit  and  that  the 
parameters  appear  to  be  proportionally  related,  we  will  next  build  a  hypothesis  function 
composed  of  two  Gaussian  functions  based  on  these  ratios.  This  function  will  have  freedom  to 
vary  in  amplitude  and  width  only.  The  widths  and  amplitudes  of  the  superimposed  Gaussians 
will  then  be  derived  from  these  meta-parameters.  In  our  construction  we  will  always  assume 
that  <jx  >  <r2  and  ax  >a2  and  we  desire  that  the  width  be  defined  as  w  =  6al  and  the  amplitude 

a  =  al+a2.  Taking  these  relations  along  with  the  ratios  /v  = ‘T’M  ancJ  Pa=aila i  we  derive 
the  Gaussian  parameters  from  w  and  a  as  follows: 


w 


(a) 


(b) 


a 

1  +  Pa 


(C) 


1  +  1/Pfl 


(d) 


(7) 


where,  again,  we  define  ^  =0.431  and  pa  =  0.551.  Our  function  then  is  constrained  by  the 
relations  in  (7)  and  takes  the  form 

-(*-*?  -o-o2 

f(x\p)  =  ale  2cr‘  +a2e  2£7:  +bx  +  c  (8) 


where  the  vector  of  free  parameters  is  now  p  =  [a  x  w  b  cf.  Because  of  the  constraints  in 

(7)  the  optimization  problem  is  once  again  well-posed  and  we  will  therefore  set  the 
regularization  parameter  C  to  zero.  After  optimizing  (1)  given  (7)  and  (8)  the  fit  functions  are 
shown  in  Figure  8.  There  it  may  be  observed  that  the  error  is  for  the  most  part  unbiased  based 
on  the  distance  from  the  center  of  the  mogul.  The  RMS  errors  for  the  right  track  are  on  the 
interval  [0.65  cm,  2.13  cm]  with  an  average  of  1.24  cm.  Likewise,  for  the  left  track  the  RMS 
error  falls  in  the  range  [0.57  cm,  2.08  cm]  with  an  average  of  1.19  cm.  Although  the 
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Figure  8.  (left)  Plot  of  the  ratio-Gaussian  best  fit  for  the  39  moguls,  (middle)  Plot  of  the 
mogul  amplitude  vs.  width  as  estimated,  (right)  Plot  of  the  error  from  the  two-Gaussian 
best  fit.  Note  scales  are  not  the  same. 

performance  is  very  slightly  worse  than  the  independent  two-Gaussian  fit  this  model  is 
appealing  because  its  shape  is  fixed  and  the  tunable  parameters  are  intuitive. 


3.6  Mogul  Model  Comparisons 

Now  that  we  have  computed  the  results  of  the  five  models,  we  desire  to  get  a  sense  as  to  how 
the  different  approaches  performed.  For  each  model,  our  optimization  produces  estimates  of  78 
moguls  (right  &  left  tracks  for  39  moguls)  for  which  we  compute  the  RMS  error.  To  gain  insight 
into  the  performance  of  each  model  we  then  compute  some  basic  statistics  of  these  78 
measures.  Table  1  shows  these  results  for  the  minimum,  maximum  and  mean  RMS  errors. 

Figure  9  shows  pictorially  how  well  our  mogul  approximations  did.  It  is  interesting  to  note  that 
regardless  of  the  chosen  model,  there  are  some  moguls  which  fit  remarkably  well  and  that  the 
worst  case  for  each  mogul  model  subjectively  looks  like  a  good  fit  as  well.  That  is  why  we  will 
now  turn  to  vehicular  response  evaluation. 


4  Vehicular  Response  Evaluation 

We  now  turn  our  attention  to  the  evaluation  of  the  modeled  mogul  shapes  from  the  point  of  view 
of  a  vehicle  dynamics  model.  In  this  section  we  will  evaluate  each  mogul  model  with  respect  to 
its  affect  on  the  response  of  an  Army  vehicle,  specifically  the  HMMWV.  For  the  original  mogul 
and  each  of  our  five  models  we  will  generate  a  short  profile  and  run  the  HMMWV  dynamics 
model  over  the  mogul  at  speeds  varying  between  8  kph  (5  mph)  and  48  kph  (30  mph)  at  8  kph 
(5  mph)  intervals  (only  four  of  these  speeds  are  reported  in  this  paper).  We  will  then  examine 
measures  from  the  dynamics  model  to  assess  the  pitch/heave  dynamics  of  the  model  under 
various  representations  of  the  same  mogul.  We  begin  with  a  description  of  the  HMMWV  model. 


Table  1 .  Statistics  of  the  RMS  error  for  39  moguls. 


Gauss 

Hann 

SVM 

2-Gauss 

Ratio-Gauss 

min  [cm] 

0.87 

1.31 

1.24 

0.57 

0.61 

max  [cm] 

2.43 

3.13 

1.64 

2.05 

2.11 

mean  [cm] 

1.62 

2.30 

1.42 

1.10 

1.21 
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Best  Median  Worst 


Figure  9.  Pictorial  representation  of  the  fit  of  the  five  different  mogul  models.  Shown  in 
the  rows  are  each  of  the  five  different  modeling  schemes  and  in  the  columns  are  the 
best,  median  and  worst  fit  for  each  mogul  model.  Axes  for  all  plots  are  [-6m,  +6m]  for 
the  horizontal  axis  and  [-0.1m,  +0.7m]  for  the  vertical  axis. 


4. 1  Dynamics  Model 

The  dynamics  model  is  a  rigid-body  model  of  the  Ml  114  High  Mobility  Multi-purpose  Wheeled 
Vehicle  (HMMWV).  Figure  10  shows  the  exterior  dimensions.  The  Gross  Vehicle  Weight 
(GVW)  of  5,970  kg  (13,160  lbs),  consisting  of  21  rigid  bodies,  includes  a  crew  weight  of  490  kg 
(1,080  lbs)  and  a  1,034  kg  (2,280  lbs)  payload.  Tire  pressures  are  set  at  276/345  kPa  (40/50 
psi)  for  the  front  and  rear  tires,  respectively. 

The  vehicle’s  tire  model  is  based  on  a  distributed  contact  method.  The  normal  deflection  and 
velocity  are  computed  based  upon  the  intersection  of  the  undeformed  tire  circle  and  terrain 
profile.  The  tire  circle  is  divided  into  vertical  ‘slices’.  Once  the  total  intersected  area  is  found,  an 
equivalent  normal  deflection  is  found  by  equating  the  intersected  area  with  that  for  the  tire  on  a 
flat  surface,  where  the  intersection  would  define  a  chord  of  the  circle.  The  point  of  application 
for  the  tire  force  is  found  through  a  weighted  average  of  the  centroids  of  the  partial  intersected 
areas  for  each  slice.  Likewise,  the  direction  of  the  force  is  found  through  a  weighted  average  of 
the  terrain  normal  vectors  associated  with  the  partial  intersected  areas  of  each  slice.  The  tire 
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•  MUDFLAPS  FROM  SIDE  VIEW  ARE  OMITTED  FOR  CLARITY 


Figure  10.  HMMWV  Exterior  Dimensions  [inches]. 

model  incorporates  vertical  force  versus  deflection  curves,  as  provided  by  the  manufacturer, 
including  the  nonlinear  effect  of  contacting  the  run-flat. 

The  suspension  model  consists  of  four  independent  double  "A"  arm  suspension  units,  one 
upper  and  lower  control  arm  for  each  wheel.  Each  "A"  arm  is  connected  to  the  chassis  by  a 
revolute  joint,  which  allows  vertical  rotation.  The  "A"  arms  are  also  connected  to  the  wheel 
spindles  by  a  spherical  joint,  which  permits  each  wheel  to  rotate  laterally  and  about  the  vertical 
axis.  This  combination  allows  each  wheel  to  move  vertically  relative  to  the  chassis  and  to  rotate 
around  a  fore-aft  or  roll  axis.  Front  tie-rods  and  rear  radius  rods,  connected  between  the 
chassis  frame  and  each  wheel  assembly,  control  the  static  toe-in  angle,  roll  steer  effect,  and  the 
wheel  steer  direction.  The  rear  radius  rods  prevent  the  rear  wheels  from  steering. 

Vertical  suspension  forces  are  provided  by  four  constant-rate  coil  springs  with  stiffness  of  350 
kN/m  (2,000  Ibs./in)  in  the  front  and  525  kN/m  (3,000  Ibs./in)  in  the  rear  and  four  telescoping 
translational  shock  absorbers  (located  inside  the  spring  elements)  represented  by  eight 
translational  springs  attached  between  the  chassis  and  the  lower  "A"  arms.  Jounce  and 
rebound  travel  is  limited  by  these  shock  absorbers.  The  shocks  provide  a  damping  force 
proportional  to  the  relative  velocity  between  the  shock  piston  rod  and  the  shock  piston  cylinder. 
The  shock  has  three  ranges  of  operation;  a  midstroke  region  and  two  hydraulic  bump  stop 
regions  (jounce  and  rebound).  The  hydraulic  bump  stops  are  built  into  the  shocks  near  the  end 
of  travel  in  each  direction  to  increase  the  damping  force  before  metal  to  metal  contact  occurs. 

The  model  was  built  in  the  Dynamic  Analysis  and  Design  System  (DADS®)  version  9.6  analysis 
tool.  There  were  six  variants  of  each  of  the  39  moguls  to  be  run  at  6  speeds  resulting  in  1 ,404 
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model  runs  in  support  of  this  analysis.  For  each  run  the  DADS  model  generated  32  channels  of 
data  reported  at  a  rate  of  100  Hz.  The  channels  reported  include  the  linear  position,  velocity, 
and  acceleration  of  the  vehicle  chassis,  acceleration  at  the  driver’s  position,  the  vehicle  pitch 
and  roll  positions  and  velocities,  the  vehicle  speed,  the  terrain  elevation  under  each  of  the  front 
wheels,  the  four  suspension  deflections,  the  four  tire  deflections  and  forces. 

4.2  Results 

In  our  analysis  we  compare  the  response  generated  by  the  mogul  model  rjt)  to  the  response 
generated  by  the  actual  mogul  profile  segment  rp(t)  by  computing  the  error  as  e(t)  =  rm(t) - rp(t) 

with  the  time  t  ranging  from  the  time  that  the  vehicle  initially  impacts  the  mogul  to  the  time 
corresponding  to  18.54  m  (50  ft.  +  130  in.)  after  impact.  We  reduce  the  error  time  history  to  a 
single  statistic  by  computing  its  RMS  value  for  selected  channels  and  only  report  relative  errors 
which  are  computed  as  erd  =RMS(e(t))/RMS(rp(t)) .  Of  the  32  channels  reported  by  the 

dynamics  model  we  will  chose  the  driver’s  seat  vertical  acceleration,  vehicle  pitch  angle,  vehicle 
pitch  angle  rate,  suspension  deflection,  tire  deflection  and  tire  force  on  which  to  compute  error 
statistics.  For  the  channels  which  are  reported  per  wheel,  we  report  the  statistic  taken  over  all 
four  corners.  These  statistics  are  shown  in  Table  2.  To  clarity  these  results  we  have  plotted  the 


Table  2.  Statistics  of  the  HMMWV  model  response  for  39  moguls. 


Speed/Channel 

Gauss 

Harm 

SVM 

2-Gauss 

Ratio-Gauss 

min 

max 

avg 

min 

max 

avg 

min 

max 

avg 

min 

max 

avg 

min 

max 

avg 

S' 

Driver's  zdd 

56.6 

94.7 

80.7 

67.9 

125.2 

113.3 

48.6 

82.0 

63.5 

23.3 

76.0 

44.8 

24.6 

139.6 

53.0 

Q- 

E 

Pitch 

5.0 

16.6 

12.3 

9.7 

23.6 

18.5 

8.0 

13.5 

10.4 

2.1 

14.3 

6.7 

2.1 

14.6 

7.5 

o 

T— 

Pitch  rate 

11.1 

32.8 

23.5 

16.6 

49.0 

34.1 

14.9 

27.0 

19.7 

5.4 

24.1 

12.9 

6.4 

25.6 

14.9 

£ 

Avg  Susp  Defl 

42.7 

53.6 

48.7 

60.2 

79.3 

70.4 

33.3 

41.0 

37.9 

23.3 

26.1 

25.0 

26.8 

31.5 

29.8 

1— 

Avg  Tire  Defl 

8.2 

9.1 

8.6 

11.0 

12.7 

11.8 

6.9 

7.1 

7.0 

4.5 

5.6 

5.0 

5.4 

6.3 

5.8 

CD 

i— 

Avg  Tire  Force 

10.4 

10.6 

10.5 

13.6 

14.8 

14.2 

8.4 

9.2 

8.8 

5.7 

7.6 

6.6 

6.7 

8.3 

7.5 

S' 

Driver's  zdd 

37.9 

87.7 

62.3 

63.7 

111.7 

97.9 

31.8 

66.8 

47.7 

17.0 

65.8 

32.2 

18.6 

96.9 

36.8 

Q. 

E 

Pitch 

4.3 

15.3 

11.1 

9.0 

23.1 

17.1 

6.9 

12.4 

9.5 

1.5 

14.6 

6.5 

2.2 

17.3 

7.3 

in 

i— 

Pitch  rate 

6.4 

29.2 

18.2 

12.4 

44.3 

27.7 

8.5 

23.0 

14.7 

3.0 

20.9 

10.2 

4.4 

21.0 

11.7 

Q. 

Avg  Susp  Defl 

37.6 

43.2 

40.1 

59.0 

64.5 

61.3 

27.5 

32.6 

29.7 

18.3 

20.2 

18.9 

21.0 

24.2 

22.3 

I- 

Avg  Tire  Defl 

12.1 

15.3 

13.8 

19.1 

22.1 

20.6 

9.3 

12.4 

10.8 

6.3 

8.1 

7.1 

7.2 

9.3 

8.3 

CM 

Avg  Tire  Force 

13.7 

17.4 

15.6 

21.8 

25.3 

23.5 

10.7 

14.2 

12.4 

7.8 

10.4 

9.1 

8.7 

11.6 

10.2 

S' 

Driver's  zdd 

21.6 

66.8 

44.7 

35.4 

82.7 

64.2 

24.9 

49.6 

36.5 

16.2 

59.8 

28.4 

18.8 

63.3 

31.3 

D. 

E 

Pitch 

5.1 

20.0 

13.2 

12.3 

30.3 

20.4 

8.2 

14.8 

11.6 

1.7 

21.5 

6.5 

1.2 

24.1 

7.8 

o 

CM 

Pitch  rate 

7.6 

28.2 

18.9 

16.6 

42.4 

29.0 

11.0 

20.2 

14.8 

2.5 

25.7 

9.0 

2.5 

29.4 

10.8 

% 

Avg  Susp  Defl 

26.5 

29.2 

27.8 

39.1 

43.9 

41.5 

20.0 

22.5 

21.2 

13.3 

15.0 

14.0 

15.5 

17.7 

16.5 

CM 

CM 

CO 

Avg  Tire  Defl 

12.1 

15.6 

13.8 

18.3 

23.0 

20.6 

9.7 

13.3 

11.5 

7.2 

8.6 

7.9 

8.0 

10.2 

9.1 

Avg  Tire  Force 

15.0 

18.4 

16.6 

21.9 

26.4 

24.0 

11.6 

15.8 

13.7 

9.4 

11.4 

10.4 

10.4 

13.1 

11.7 

S' 

Driver's  zdd 

22.9 

61.3 

38.1 

18.5 

75.1 

52.8 

19.9 

64.8 

40.0 

9.2 

68.9 

26.9 

9.5 

74.7 

28.2 

Q. 

E 

Pitch 

3.9 

21.0 

10.8 

4.8 

27.8 

16.9 

3.3 

21.8 

12.0 

1.0 

29.7 

7.6 

1.7 

32.0 

7.7 

in 

CM^ 

Pitch  rate 

5.8 

26.0 

15.1 

6.8 

38.4 

23.5 

5.2 

29.0 

16.1 

1.6 

36.6 

9.7 

2.5 

39.3 

10.2 
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Figure  11.  Plots  showing  the  relative  errors  vs.  speed  for  each  of  the  five  mogul 
models.  The  plots  show  the  driver’s  vertical  acceleration,  pitch,  pitch  rate,  suspension 

deflection,  tire  deflection  and  tire  force  respectively.  The  a,  v,  and  o  symbols 


represent  the  minimum,  maximum  and  mean  respectively.  The  x  symbol  shows 

individual  relative  errors  for  each  of  the  moguls.  The  bottom  three  plots  are  shown  for 
each  of  the  four  corners.  For  each  speed  the  mogul  models  from  left  to  right  are  the 
Gaussian,  Planning  Window,  SVM,  2-Gaussian  and  ratio-Gaussian.  Note  that  plots  are 
at  different  scales. 


same  statistics  in  Figure  11.  There  we  observe  that  typically  the  best  performer  is  the  two- 
Gaussian  model  followed  closely  by  the  ratio-Gaussian  model,  the  next  best  model  appears  to 
be  the  SVM  followed  by  the  Gaussian  model.  The  worst  performer  is  the  Planning  window.  We 
may  also  examine  Figure  1 1  from  the  perspective  of  speed.  It  appears  that  for  the  vehicle 
sprung  mass  there  appears  to  be  a  negligible  to  slight  improvement  in  the  relative  error  with 
increasing  speed.  On  the  other  hand  for  the  tire  deflection  and  force,  there  is  a  slight  increasing 
trend  with  vehicle  speed.  The  worst  performing  channel  is  the  driver’s  vertical  acceleration.  We 
suggest  that  this  is  the  case  because  of  the  roughness  in  the  original  profile  which  is  not  present 
in  the  mogul  model.  Acceleration  signals  are  particularly  sensitive  to  the  high-frequency  nature 
of  the  roughness  and  thus  the  relative  roughness/smoothness  is  highlighted  in  acceleration 
signals. 

Since  the  two-Gaussian  and  ratio-Gaussian  models  perform  similarly  it  is  reasonable  to  ask: 
which  is  the  better  model?  Although  the  two-Gaussian  function  performs  slightly  better,  it  does 
not  define  a  shape  for  the  prototypical  mogul  due  to  its  mathematical  structure.  Our  desire  is  to 
obtain  a  prototypical  shape  which  may  be  parameterized  by  width  and  amplitude.  In  this  case 
we  choose  the  ratio-Gauss  function  to  be  the  best  representation  of  the  39  moguls  that  we  have 
examined.  Given  this  conclusion,  we  illustrate  (in  Figure  12)  the  performance  of  this  model  for 
the  median  performer,  mogul  number  39.  There  we  observe  time  histories  of  the  channels  used 
in  this  study.  By  observing  the  driver’s  vertical  acceleration,  it  is  evident  that  the  high-frequency 
components  of  the  signal  are  largely  contributing  to  the  error.  Some  of  this  high  frequency 
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Figure  12.  Plots  showing  the  reference  (blue)  and  the  model  (green)  responses  for  the 
ratio-Gauss  mogul  at  32.2  kph  (20  mph),  shown  are  responses  from  mogul  39  which 
represents  the  responses  whose  error  are  closest  to  the  median  error. 


behavior  may  be  observed  in  the  tire  deflection  and  forces  as  well.  We  should  note  that  at 
higher  speeds,  the  shape  of  the  mogul  dominates  the  vehicle  response,  while  at  the  lower 
speeds,  the  surface  roughness  tends  to  dominate  the  response.  Although  the  responses  shown 
in  Figure  12  subjectively  appear  to  be  very  good,  they  would  be  worse  at  the  slower  speeds. 


5  Conclusion 

In  this  paper  we  have  discussed  modeling  the  shape  of  discrete  features  on  an  Army  endurance 
test  course.  We  illustrated  why  existing  statistical  techniques  are  not  well  suited  for 
characterizing  and  synthesizing  these  features.  We  offer  five  functions  as  candidates  with 
which  to  model  the  shape  of  the  moguls.  These  shape  functions  varied  in  the  number  of  free 
parameters  available.  The  Gaussian  and  the  Planning  window  shapes  were  parameterized  by 
width  and  amplitude,  the  SVM  function  technically  has  161  free  parameters,  although  most 
vanish  in  the  optimization.  The  two-Gaussian  function  has  four  free  parameters  and  the  ratio- 
Gauss  has  two  free  parameters.  We  evaluated  the  quality  of  each  model  by  computing  the 
RMS  elevation  error  of  each  fit  and  by  examining  the  dynamic  response  of  a  PIMMWV  model. 
For  the  PIMMWV  model  we  examined  six  different  types  of  signals  related  to  the  chassis, 
suspension  and  tires.  We  evaluated  the  quality  of  the  mogul  model  by  again  computing  the 
RMS  error  over  each  reference  signal  for  four  different  speeds.  After  examining  the  response 
errors,  we  concluded  that  the  two-Gaussian  and  ratio-Gaussian  functions  give  the  best 
approximation  of  the  moguls.  Among  these  two  functions  we  conclude  that  the  ratio-Gauss 
function  is  preferred  because  it  has  fewer  free  parameters  and  has  a  definite  shape.  We 
therefore  conclude  that  the  39  moguls  on  this  test  course  are  best  modeled  by  (7)  and  (8)  with 
pa  =0.431,  pa=  0.551  and  the  distribution  of  width  and  amplitude  given  by  the  distribution 
shown  in  Figure  8. 
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The  analysis  in  this  paper  was  conducted  for  the  profile  as  measured  in  2002.  This  profile  has 
been  periodically  surveyed  by  the  Army’s  profilometer  for  a  decade.  Future  work  may  include 
an  analysis  of  how  these  moguls  change  over  time.  Also,  there  are  other  potential  candidate 
functions  which  were  not  evaluated  in  this  study.  Future  work  may  include  the  evaluation  of 

functions  such  as  sech(-)  and  ^tanh(.*)  =  l-tanh2(;c) . 
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